/**********************************************************************
Compile : gcc -O3 -o generate3 check.c generate3.c
Usage : ./generate3 mode mode1 mode2 pnflag N L M n seed       
mode  : 0 common; 1 cycled; 2: triggled cycle
mode1 : 0 disallow self-activation; 1 allow self-activation
mode2 : Additional Requirements: 0 no requirements; 1 restriction I; 2 restriction II; 3 restriction III
pnflag: 1 print the candidate processes
N     : number of the nodes                
L     : number of the steps        
M     : the number of "1" in the process  
n     : the number of the attempts         
seed  : random seed
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

void print_help ()
{
	printf ("Usage         : ./generate3 mode mode1 mode2 pnflag N L M n seed\n");
	printf ("                (./generate3 0 0 0 1 11 12 41 10000 $RANDOM)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("pnflag        : 1 print the candidate processes\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n");
	printf ("M             : the number of \"1\" in the process\n");
	printf ("n             : the number of the attempts\n");
	printf ("seed          : random seed\n\n");
	printf ("Note          : This program is special optimized for sampling under restriction I\n");
}

int bricknum;

int check_bricks (int L, int brick[L], int mode)
{
	int i;

	if (mode != 2)
	{
		if (brick[0] == 0 && brick[1] == 1 && brick[2] == 0)
			return 0;
		if (brick[0] == 1 && brick[1] == 0 && brick[2] == 1)
			return 0;
	}

	for (i = 1; i < L - 3; i++)
	{
		if (brick[i] == 0 && brick[i + 1] == 1 && brick[i + 2] == 0)
			return 0;
		if (brick[i] == 1 && brick[i + 1] == 0 && brick[i + 2] == 1)
			return 0;
	}

	if (mode == 1)
	{
		if (brick[L - 3] == 0 && brick[L - 2] == 1 && brick[0] == 0)
			return 0;
		if (brick[L - 3] == 1 && brick[L - 2] == 0 && brick[0] == 1)
			return 0;
		if (brick[L - 2] == 0 && brick[L - 1] == 1 && brick[1] == 0)
			return 0;
		if (brick[L - 2] == 1 && brick[L - 1] == 0 && brick[1] == 1)
			return 0;
	}

	return 1;
}

void combination1 (int L, int n, int brick[L], int bricks[][L], int bricks_np[], int mode)
{
	int i;

	if (n == L - 1)
	{
		if (mode == 1)
			brick[L - 1] = brick[0];
		else
			brick[L - 1] = brick[L - 2];

		if (mode == 2)
			brick[0] = brick[L - 1];

		if (check_bricks (L, brick, mode) == 1)
		{
			memcpy (bricks[bricknum], brick, sizeof (int) * L);
			if (mode == 2)
				L --;
			bricks_np[bricknum] = 0;
			for (i = 0; i < L - 1; i++)
			{
				bricks_np[bricknum] += brick[i];
			}
			if (bricks_np[bricknum] > 0 && bricks_np[bricknum] < L - 1)
				bricknum ++;
		}
	}
	else
	{
		brick[n] = 0;
		combination1 (L, n + 1, brick, bricks, bricks_np, mode);
		brick[n] = 1;
		combination1 (L, n + 1, brick, bricks, bricks_np, mode);
	}
}

void combine_dist (int N1, int N2, int dist[N2 + 1], double Dist[N1][N1 * N2 + 1], int i)
{
	int j, k;

	if (i == 0)
	{
		for (j = 0; j <= N1 * N2; j++)
		{
			Dist[i][j] = 0.0;
		}
		for (j = 0; j <= N2; j++)
		{
			Dist[i][j] = dist[j];
		}
	}
	else
	{
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist[i][j] = 0.0;
			if (j <= N2)
			{
				for (k = 0; k <= j; k++)
				{
					Dist[i][j] += dist[k] * Dist[i - 1][j - k];
				}
			}
			else
			{
				for (k = 0; k <= N2; k++)
				{
					Dist[i][j] += dist[k] * Dist[i - 1][j - k];
				}
			}
		}
	}
}

void combine_dist2 (int N1, int N2, int dist[N2 + 1], double Dist[N1 * N2 + 1], int i)
{
	double Dist2[N1 * N2 + 1];
	int j, k;

	if (i == 0)
	{
		for (j = 0; j <= N1 * N2; j++)
		{
			Dist[j] = 0.0;
		}
		for (j = 0; j <= N2; j++)
		{
			Dist[j] = dist[j];
		}
	}
	else
	{
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist2[j] = 0.0;
			if (j <= N2)
			{
				for (k = 0; k <= j; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
			else
			{
				for (k = 0; k <= N2; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
		}
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist[j] = Dist2[j];
		}
	}
}

void setup_count (int N, int L, double Count[N * (L - 1) + 1])
{
	int i, j, count[L];

	count[0] = count[1] = 1;

	for (i = 1; i < L - 1; i++)
	{
		count[i + 1] = 1;
		for (j = i; j > 0; j--)
		{
			count[j] = count[j] + count[j - 1];
		}
	}

	count[0] = count[L - 1] = 0;

	for (i = 0; i < N; i++)
	{
		combine_dist2 (N, L - 1, count, Count, i);
	}
}

void setup_possibility (int N, int L, int n, int m, int bricks_np[], int count[], double Count[N][N * (L - 1) + 1], double possibility[])
{
	int i;
	double p[L];

	if (n == N - 1)
	{
		for (i = 0; i < L; i++)
		{
			p[i] = 0.0;
		}
		p[m] = 1.0;
	}
	else
	{
		for (i = 0; i < L; i++)
		{
			p[i] = count[i] * Count[N - 2 - n][m - i] / Count[N - 1 - n][m];
		}
	}

	for (i = 0; i < bricknum; i++)
	{
		possibility[i] = p[bricks_np[i]] / count[bricks_np[i]];
	}

	for (i = 1; i < bricknum; i++)
	{
		possibility[i] += possibility[i - 1];
	}

	/*
	printf ("Round : %3d %d\n", n, m);
	for (i = 0; i < bricknum; i++)
	{
		printf ("%3d : %3d %lf\n", i, bricks[i][0], possibility[i]);
	}
	*/
}

int roulette (double possibility[], int low, int upper, double p)
{
	if (upper == low + 1)
		return upper;
	else
	{
		int med;

		med = (low + upper) / 2;

		if (p > possibility[med])
			return roulette (possibility, med, upper, p);
		else
			return roulette (possibility, low, med, p);
	}
}

void generate_pathway (int N, int L, int M, int pathway[L][N], int bricks[][L], int bricks_np[], int count[], double Count[N][N * (L - 1) + 1], int mode)
{
	int n, i, j, m, select;
	double p, possibility[bricknum];

	m = M;
	for (i = 0; i < N; i++)
	{
		if (mode == 2 && i == 0)
		{
			pathway[0][0] = 1;
			for (j = 1; j < L; j++)
			{
				pathway[j][0] = 0;
			}
			continue;
		}

		if (mode == 2)
			setup_possibility (N - 1, L - 1, i - 1, m, bricks_np, count, Count, possibility);
		else
			setup_possibility (N, L, i, m, bricks_np, count, Count, possibility);

		p = drand48 ();
		if (p < possibility[0])
			select = 0;
		else
			select = roulette (possibility, 0, bricknum, p);

		m = m - bricks_np[select];

		for (j = 0; j < L; j++)
		{
			pathway[j][i] = bricks[select][j];
		}
	}
}

int check_duplicate (int N, int L, int pathway[L][N])
{
	int i, j, flag;

	for (i = 0; i < L - 2; i++)
	{
		flag = 1;
		for (j = 0; j < N; j++)
		{
			if (pathway[i][j] != pathway[i + 1][j])
			{
				flag = 0;
				break;
			}
		}
		if (flag == 1) break;
	}

	return flag;
}

double factorial (int n)
{
	double F;

	if (n == 0)
		return 0;

	F = n * (log (n) - 1) + 0.5 * log (2.0 * M_PI * n) + 1.0 / (12.0 * n) - 1.0 / (360.0 * n * n * n);

	return F;
}

int main (int argc, char *argv[])
{
	int i, j, k, m, flag, pnflag;
	int N, L, M, n, mode, mode1, mode2;
	long int seed;

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	mode = atoi (argv[1]);
	mode1 = atoi (argv[2]);
	mode2 = atoi (argv[3]);
	pnflag = atoi (argv[4]);
	N = atoi (argv[5]);
	L = atoi (argv[6]);
	M = atoi (argv[7]);
	n = atoi (argv[8]);
	seed = atoi (argv[9]);
	srand48 (seed);
	printf ("Random Seed : %ld\tAttempts : %d\n", seed, n);

	int pathway[L][N];

	/* Setup bricks */
	bricknum = 1;
	for (i = 0; i < L - 1; i++)
	{
		bricknum *= 2;
	}

	if (mode == 2)
	{
		bricknum /= 2;
	}

	int brick[L], bricks[bricknum][L], bricks_np[bricknum];

	bricknum = 0;
	if (mode == 2)
		combination1 (L, 1, brick, bricks, bricks_np, mode);
	else
		combination1 (L, 0, brick, bricks, bricks_np, mode);
	/* End of Setup bricks */

	if (mode == 2)
	{
		N --; L --;
	}

	/* Setup the number of all possible processes for each given activity */
	/* count[i]: the number of bricks with activity = i */
	int count[L];

	for (i = 0; i < L; i++)
		count[i] = 0;

	for (i = 0; i < bricknum; i++)
		count[bricks_np[i]] ++;

	/* Count[i][j]: the number of none quick-flip ixL processes with activity = j */
	double Count[N][N * (L - 1) + 1];

	for (i = 0; i < N; i++)
	{
		combine_dist (N, L - 1, count, Count, i);
	}

	/* Count2[i]: the number of processes with activity = i (all-zero node is excluded) */
	double Count2[N * (L - 1) + 1];

	setup_count (N, L, Count2);
	/* End of Setup */

	m = 0;

	if (mode == 2)
	{
		N ++ ; L ++;
	}

	for (i = 0; i < n; i++)
	{
		generate_pathway (N, L, M, pathway, bricks, bricks_np, count, Count, mode);

		if (check_duplicate (N, L, pathway) == 1)
			continue;

		flag = network_inferred (N, L, pathway, mode, mode1, mode2);

		if (flag == 1)
		{
			if (pnflag == 1)
			{
				printf ("c candidate process %d\n", M);
                                
				for (j = 0; j < L; j++)
				{
					for (k = 0; k < N; k++)
					{
						printf ("%d ", pathway[j][k]);
					}
					printf ("\n");
				}
			}

			m ++;
		}
	}

	if (mode == 2)
		N --;

	printf ("Total %d candidate processes with activity %3d are obtained.\n", m, M);
	printf ("Estimated feasible processes with activity %3d : %lg\n", M, Count[N - 1][M] * m / n);
	printf ("Estimated possibility : %lg\n", (double) (Count[N - 1][M] * m) / (double) (Count2[M] * n));

	exit (0);
}
